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We examine how the properties of the Kondo insulators change when the symmetry of the underlying crys- 
tal field multiplets is taken into account. We employ the Anderson lattice model and consider its low-energy 
physics. We show that in a large class of crystal field configurations, Kondo insulators can develop a topolog- 
ical non-trivial ground-state. Such topological Kondo insulators are adiabatically connected to non-interacting 
insulators with unphysically large spin-orbit coupling, and as such may be regarded as interaction-driven topo- 
logical insulators. We analyze the entanglement entropy of the Anderson lattice model of Kondo insulators by 
evaluating its entanglement spectrum. Our results for the entanglement spectrum are consistent with the surface 
state calculations. Lastly, we discuss the construction of the maximally localized Wannier wave functions for 
generic Kondo insulators. 



PACS numbers: 71.27.+a, 75.20.Hr, 74.50.-M 

I. INTRODUCTION 

Topological insulators are a novel class of materials in 
which strong spin-orbit interaction leads to the inversion of the 
band gap (See Ref.Q]and|2]and references therein). In 3D, this 
inversion results in chiral metallic surface states due to a for- 
mation of a single Dirac cone inside the gap 3 -"^. Among mate- 
rials which exhibit this behavior, for example, HgTe, Bi2Se3, 
Bi 1 _ 2 .Sb :I ., Bi 2 Te 3 and TlBiTe 2 , the chiral structure of the sur- 
face states has been confirmed experimentally^^] 

The emergence of surface modes in topological insulators 
is a band structure effect which can be understood without in- 
voking interactions. There is great current interest in the pos- 
sibility of interaction driven topological phenomena. Up until 
now there are no experimental examples of interaction-driven 
topological insulators that preserve time-reversal symmetry. 
However, several theoretical proposals have been put forward: 
2D topological insulators via spontaneous symmetry break- 
ing in bilayer graphene and optical lattice systems^EIl topo- 
logical Mott insulating phase in Ir-based pyrochlore oxides 
A 2 Ir 2 07 with A=Nd,P r 6 " 19 [ Kondo insulators with the most 
salient example of SmBg 20 and insulating behavior in filled 
skutteridites 21 . In this paper we focus on general principles 
governing the emergence of chiral metallic states in Kondo 
insulators. Throughout this paper, we use the term "Kondo 
insulator" in its broadest sense, including both mixed valent 
materials 22 such as SmBg and YbBi 2 and those in the more 
localized limit, such as Ce3Bi4Pt3. 

Kondo insulators are a type of heavy fermion material, first 
discovered forty years age 23 , in which highly renormalized /- 
electrons hybridize with conduction electrons to form a com- 
pletely filled band of quasiparticles with excitation gaps in the 
millivolt range^H^Zl. Because Kondo insulators appear as a 
result of strong interactions, one might think that their ex- 
citations and their ground-states are adiabatically connected 
to trivial non-interacting band insulators 28 . However, before 
jumping to this conclusion one needs to be careful, for in the 



renormalization process the width of the heavy electron bands 
drops far below the characteristic size of the spin-orbit inter- 
actions, driving the physics to a new fixed point characterized 
by infinite spin-orbit coupling in the localized bands. Indeed, 
we shall show that topological Kondo insulators are adiabati- 
cally connected to non-interacting topological insulators with 
an unphysically large value of the spin-orbit coupling, and in 
this sense, they are interaction-drive insulators. 

One of the most important features of the /-electron sys- 
tems in general and Kondo insulators in particular is that the 
/-electron states are classified with respect to their momen- 
tum k, total angular momentum J and its z-axis component 
M, while conduction electron states are described by a mo- 
mentum and a spin a. When an f-electron escapes into the 
conduction sea, it hybridizes with a spin-orbit coupled Wan- 
nier state of the conduction electrons that has the same sym- 
metry as the f-state. The spin-orbit coupled Wannier states 
of the conduction electrons are then decomposed in terms of 
plane-wave states and this gives rise to momentum-dependent 
form factors with symmetries that are uniquely determined by 
the local symmetry of the f-states. In this way, the form- 
factors encode the effect of the strong spin-orbit coupling. 
More importantly, these form factors also define the underly- 
ing symmetry of the hybridization amplitude and gap which is 
develops below the "Kondo temperature" Tk at which heavy 
quasi-particles develop. One of the key properties of the spin- 
orbit coupled f-state, is an odd-parity wavefunction. It is the 
protected odd-parity of the f-states that provides the driving 
force for the formation of topological insulating states. 

The dimension of the form factor matrix is determined by 
the degeneracy of the underlying ground state /-ion multiplet. 
In the crystalline environment the (2 J + 1) multiplet degener- 
acy is lifted by the crystalline fields. For half-integer values of 
J the lowest possible degenerate multiplet is a Kramers dou- 
blet, which means that the form factor is a two-dimensional 
matrix. For the integer values of J the crystal field can fully 
lift degeneracy of the multiplet. This situation corresponds to 



2 



the non-magnetic state of the /-ion and currently there are no 
known examples of such Kondo insulators. Thus, in this arti- 
cle we will only focus on the magnetic ions with half-integer 
values of the total angular momentum. 

As mentioned above, the symmetry of the lowest lying mul- 
tiplet determines the symmetry of the hybridization ampli- 
tude. Generically, two possible scenarios can arise depend- 
ing on whether the hybridization contains nodes or not. For 
a small, but important subset of these systems, the hybridiza- 
tion contains nodes. In this case, the Kondo insulating state 
is replaced by a heavy semi-metal with a pseudogap, as in the 
case of CeNiSn or CeRhSb. If the nodes correspond to touch- 
ing of the two non-degenerate bands with linear dispersion, 
the system becomes a Weyl semi-metal, where topologic ally- 
protected surface modes emerge 2 ^. Note, that the lifting of 
degeneracy can only happen due to onset of magnetic order. 
The magnetic moments may appear as a result of incomplete 
screening similarly to what happens in CeCoIn5, for example. 

In our previous work, we demonstrated within a mean-field 
model, that in the large class of systems without nodes in the 
hybridization gap, Kondo insulators can develop topological 
insulating ground-states®. In this paper we develop this idea 
in detail, providing mathematical details of the construction of 
the wavefunction and explicitly computing the surface modes 
for a topological Kondo insulator. 

Although attempts to establish the general principle deter- 
mining the relative position of the crystal field multiplets have 
been mad d 30 " 32 ! but no such principle has yet been discovered. 
Experimentally, however, the symmetry of the lowest lying as 
well as excited multiplets can be detected, for example, by 
inelastic neutron scattering spectroscopy 33 . Nevertheless, by 
assuming a specific symmetry of the ground state multiplet 
one is able to theoretically predict the physical properties of a 
Kondo insulator. In addition, it provides multiple ways to ver- 
ify them experimentally. To be more precise, the presence of 
the chiral states on the surface of Kondo insulators will allow 
one to indicates unambiguously the symmetry of the lowest 
lying multiplet. 

Apart from going beyond the brief description of topolog- 
ical Kondo insulators reported in Ref. |20] we will discuss 
the topological properties of the eigenfunctions of the model 
Hamiltonian describing the Kondo insulators. We will start 
with a short review of the model and recently obtained results. 
We evaluate the entanglement spectrum for simplest model of 
the Kondo insulator corresponding to the nearest neighbors 
tight-binding approximation for the conduction bands and a 
Kramers doublet. We also discuss the choice of the proper ba- 
sis for construction of the maximally localized Wannier func- 
tions for the Kondo insulators. We show that Wannier func- 
tions can be constructed on the basis composed of the linear 
combination between the conduction and /-electron states. 
Finally, we provide a short review of available experimental 
data which points towards the existence of the chiral surface 
states in Kondo insulators. 



II. ANDERSON LATTICE MODEL 

We begin with writing down the model Hamiltonian to de- 
scribe the physics of the Kondo insulators. In what follows 
we will consider the most general case by assuming that there 
are N c conduction bands, so that the Hamiltonian describing 
conduction electrons is 



H, 



= ££fc*£M? 



(i) 



1 = 1 k,tr 



where £jk is the dispersion of the kh band of conduction elec- 
trons, a is the spin index and c^J is a conduction electron 
creation operator. Consequently, the Hamiltonian which de- 
scribes the /-electrons is: 



JYr 



H f ~ EE e /r/]a/j"a + U ^ fiafiafia'fi, 



(2) 



where fj a creates an /-electron on site j in a state a of a low- 
est lying multiplet Np -degenerate multiplet denoted by T (see 
below), Cf is the /-electron energy and U > is the strength 
of the Hubbard interaction between the /-electrons. We em- 
phasize that index a is not a spin index due to the presence of 
the strong spin-orbit coupling. Generally states belonging to 
the multiplet T are described by the total angular momentum 
J and z-component M or some linear superposition of those 
states and in the second term of Eq. d2l, the summation is 
restricted to a ^ a'. 

Finally the term describing how electrons in N c conduction 
bands are hybridized with localized /-electrons is 



N c N r 



1=1 j,a=l 



(3) 



Here V^j a is a non-local hybridization matrix element be- 
tween the conduction electrons in kh band and localized /- 
electrons. Thus, the periodic Anderson model Hamiltonian, 
which is the basis for our subsequent discussion, reads: 



Hp am — H c + Hf + Hfr 

r(0 



(4) 



The hybridization matrix elements V ity j a can be written as 
follows : 



(5) 



ker 



where Vi is the hybridization amplitude and the form factors 
[^rk]ao- are (2 J + 1) x 2 dimensional matrices given by: 



[$rk]a ff = (kTa\ka) 



(6) 



Since in this paper we will be discussing the materials when 
an /-ion is in the valence state with J = 5/2 (f 1 for cerium 
or / 3 for samarium) it follows 
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3m,-a)F^_ CT (k) (7) 
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and 



Y 3 



(k) 



Z 



Y 3 



(R)e 



jk-R 



(8) 



R^O 



is a tight-binding generalization of the spherical Harmonics 
that preserves the translational symmetry of the hybridization, 
$(k) = <£>(k + G), where G is reciprocal lattice vector. Here, 
R are the positions of the Z nearest neighbor sites around the 
magnetic ion. Note, that deriving ([6]) we have assumed that 
the symmetry of the conduction electron amplitude coincides 
with the symmetry of the /-ion multiple ! 34 ! 35 !. Consequently, 
we treat the system with only one hybridization channel. Now 
let us recall the definition of the form factors: 



5/2 

[$kW= £ (ka\JM)(JM\ka) 

M=-5/2 



(9) 



where ( JM \ka) is a (2 J + 1) x 2 matrix whose elements are 

given by a\j 2 !f ^-Y 3 { _ La (\i). The elements of the matrix 

(ka\ JM) are determined by the specific choice of the /-ion 
multiplet and the corresponding wave-functions denoted by 
\Ta). As we have already mentioned above, we will focus our 
discussion on the case of /-ion with J = 5/2. This situation 
is relevant for all known /-electron Kondo insulators. Conse- 
quently, in a cubic crystal field environment, the magnetic ion 
multiplet is split into a doublet 



\r[ c) ±) 



± 



±- 



(10) 



and a quartet 



where the mixing angle /3 defines orientation of the cor- 
responding states. In an orthorhombic environment, the 
Kramer's doublets are generally described by a linear super- 
position of all three wave- functionPEH 

\T(° rth °)±) = u \± 1/2) + v\ T 3/2) +w\± 5/2). (15) 

Having provided the scheme for the computation of the form- 
factors we proceed with the discussion of the low-energy 
properties of our model Q of Kondo insulators. 



III. LOW-ENERGY THEORY FOR Ce-BASED KONDO 
INSULATORS 

The low-energy properties of the model Q are described 
in terms of renormalized quasiparticles formed via strong hy- 
bridization between the c— and /— states and on-site repul- 
sion U between the /-electrons. In the regime where the / 
states are predominantly localized, U ~ W (W is the band- 
width), we can neglect the momentum dependence of the /- 
electron self-energy £/(k, u) ~ 

Below we discuss the topological properties of the effective 
low-energy model. To make our discussion more tractable, we 
will consider separately several experimentally relevant cases. 
In what follows we discuss the simplest case of the single con- 
duction band and Kramers doublet as a ground state multiplet 
of the magnetic ion. This is done with an eye toward the trans- 
port experiments on the Ce-based Kondo insulator s 36 ! 37 ! 



A. single conduction band hybridized with the Kramers 
doublet: Ce-based Kondo insulators 
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so that for this case the matrix (ka\ JM) is 
(kai\JM) = 







(11) 



(12) 



(13) 



In a tetragonal crystal field environment, relevant for Ce- 
based Kondo insulators, the Ce multiplet is split into three 
doublets: 



|rf±) 
|rf ±) 
|if ±) 



I ±1/2), 

cos(^) I T 3/2) + sin(/3) | ± 5/2) , (14) 
sin(/3)|T3/2)-cos(/3)|±5/2), 



In order to derive an effective low-energy model for Kondo 
insulators, we first introduce the following correlation func- 
tions for c- and /-electrons: 



G cc (k,r) = -<f T {c ka (r)c t k(7 (0)}), 
G / /(k,r) = -(f T {/ ka (r)/t a (0)}}, 



(16) 



By writing down equations of motion for the c-operators with 
the Hamiltonian Q and going into Matsubara frequency rep- 
resentation we derive the following relation: 



G cc (k,iw) = Gi° c ) (k, iW ) 



\V\ 2 Al 
(iuj - £k) : 



•G//(k,i 



(17) 



with A£ = |Tr[$^ k <l> rk ] and Gfc (k, iu) is a conduction 
electron propagator in the absence of interactions. If we de- 
note the /-electron self-energy by E/ (k, oj) and keep in mind 
that this self-energy appears as a result of Hubbard correla- 
tions only, then it follows: 



(18) 



Next we assume that the self-energy is very weakly dependent 
on momentum, S^(k, iuj) ~ Y.f(kp,iuj) (kp is the conduc- 
tion electron's Fermi momentum) and expand it to the lowest 



order in Matsubara frequency: 
S/(k, ioS) ~ Sf (fop, iw) + jo; 



d(iuj) 



(19) 



Taking into account expressions ( 18|19) , for the correlators 
we find 



G cc (k, iw) 



Gf f (k,iw) = 



(iu,-&)(iu,-e f )-\V\*Al' 

( lw -e k )^- £/ )-|vi 2 A2' 



(20) 



where £k = — 2t X) a =a; y z cos ^ a * s tne ^ >are s P ect; nim of 
conduction electrons taken relative to the chemical potential, 
ef = Z[tf + S/(0)] is the renormalized /-level, V = \fZV 
and Z = (1 — dY,f(kF,u))/dw)~= . These propagators cor- 
respond to the following effective HamiltoniarPJ 



rk 



Ik 



(21) 



Here 1 denotes the unit 2x2 matrix. The KI is formed 
if the chemical potential of the quasiparticles lies inside the 
hybridization gap, separating the two bands with the spectra 



E±(k) = i[a+e/±y(a-e/) 2 + 4 va, 

To discuss the topological properties of our effective model 
for the KI (21) , we need to consider separately the form fac- 
tors for different Vs. It is convenient to distinguish these 
states according to their orbital symmetry parameterized by 
the index a = 1,2,3 and the pseudo-spin quantum number 
(a = ±jSl. Hence, we have f$ ± \0) = \ ± 1/2), / 2 f ± |0) = 

|±3/2),and/] ± |0) =|±5/2). 

The momentum-dependence of 
A a (k) follows from Eq. Q. 

Ai(k) = [12cos(2(9) + 5(3 + cos(460)] 1/2 , A 2 (k) = 



the hybridization gap 
At small momenta k, 



sin6>|[17+15cos(260] 1/2 ,and A 3 (k) 



£ sin 2 , 



where 6 and 4> define the direction of the unit vector k, as- 
sociated with the point on the Fermi surface. Note that the 
hybridization gap has a line of nodes along the z-axis for the 
shapes a = 2,3, but generic combinations of all three form- 
factors characteristic of contain no nodes. The key results 
of this Section are most simply illustrated using the nodeless 
a = 1 Kramers doublet as the ground-state of the magnetic 
ion. 

To analyze the topology of the bands we use the fact that 
topology is invariant under any adiabatic deformation of the 
Hamiltonian. We begin our study with a tight-binding model 
for a KI on a simple cubic lattice. Our choice of hybridiza- 
tion ensures that the mean-field Hamiltonian (Eq. |2T) is a 
periodic function satisfying "H e //(k) = H e ff(k + G). The 
technical analysis is readily generalized to more complicated 
cases as discussed below. The most important element of 
the analysis is the odd parity form factor of the / electrons, 
$ a (k) = — <I> a (— k). This parity property is the only essential 
input as far as the topological structure is concerned. 



(a) 



(1;000) 



(b) 



(0;111) 



FIG. 1: Two topological classes can be realized in our model of 
Kondo insulators for ej < 2t (see text). The first class with in- 
dex v = (1; 000) corresponds to a strong topological insulator and is 
realized when £/ < — 2t. The second class with index v = (0; 111) 
is realized for —2t < e < 2t. When the renormalized position of the 
/-level is at the boundaries, £/ = ±2£, ±6t the system is metallic. 



B. calculation of topological indices 

In Ref. [38] Fu and Kane demonstrate that in an insulator 
with time-reversal and space-inversion symmetry , the topo- 
logical structure is determined by parity properties at the 
eight high-symmetry points, k m , in the 3D BZ which are 
invariant under time-reversal, up to a reciprocal lattice vec- 
tor: k^ = — k* n + G. In our case, these symmetries re- 
quire that Ueff(k) = PHeffi-kjP- 1 and"H e// (k) T = 
TH e ff(— k)T _1 , where the parity matrix P and the unitary 
part of the time-reversal operator T are given by 



P 



1 



-1 



T 



f i&2 



(22) 



where a 2 is the second Pauli matrix. For any space-inversion- 
odd form factor, it follows immediately that $ a (k) = at 
a high-symmetry point. Hence, the Hamiltonian at this high 
symmetry point is simply H e ff(k m ) = (£ k ^ + e/)//2 + 
(£k* ^ £ /)-P/2> where I is the four-dimensional identity ma- 
trix. 

The parity at a high symmetry point is thus determined 
by <^m — sgnfec* — Four independent Z 2 topologi- 
cal indices (i>q\ v\, 1/2, ^3) 39 , one strong (a = 0) and three 
weak indices (a = 1,2,3) can be constructed from 5 m : 
(i) The strong topological index is the product of all eight 

S m 's: 7 STI = (-1)*° = n 8 m = ±1; (ii) by set- 

ting kj — (where j — x,y,and z), three high-symmetry 
planes, Pj — {k : kj — 0}, are formed that contain four high- 
symmetry points each. The product of the parities at these 
four points defines the corresponding weak-topological index, 
= (-1)"« = II S m = ±1, a = 1,2,3 with inte- 

gers corresponding to the axes x,y and z. The existence of 
the three weak topological indices in 3D is related to a Z 2 
topological index for 2D systems (a weak 3D TI is similar to 
a stack of 2D Z 2 topological insulators). Because there are 
three independent ways to stack 2D layers to form a 3D sys- 
tem, the number of independent weak topological indices is 



ja 

i WTI 
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C. surface state calculation 

In addition to the method discussed by us in RefP^, this can 
be proven by the direct calculation of the band spectrum to- 
gether with the calculation of the entanglement entropy (see 
below). In what follows, we will assume that the /-electrons 
have very weak hole-like dispersion, i.e. Ef — > £/k = 
2tf coski + Hf, where t/ = O.lt and fit is a chemical 

i— x,y ,z 

potential. This implies, in particular, the the boundary sepa- 
rating the WTI and STI are now given by /.i' c = ±2(t + tf) 
Here, as before, the value of fi' c is taken relative to the chemi- 
cal potential of the conduction electrons. 

In order to demonstrate that there is a metallic surface state 
in the spectrum described by the Hamiltonian (21 1 we con- 



sider a stack of N = 30 planes along the z-direction and di- 
agonalize the Hamiltonian. The resulting Hamiltonian matrix 
has blocks along the diagonal, which describe the hopping and 
hybridization within each plane and the off-diagonal parts de- 
scribing the hopping and hybridization between the planes. I 
For the set of the parameters corresponding to the strong topo- 
logical Kondo insulator we compute the spectrum numerically 
and show the results on Fig. [2] For simplicity we have chosen 
the model form factor, given by: 



y(sin k x a x + sin k y a z ), within the planes, 
$ = { iV z a z , between the planes (upwards), 

-iV z a z , between the planes (downwards). 



(23) 



FIG. 2: Single particle band spectrum governed by the mean field 
Hamiltonian \H\ along the a>axis, k — k x ,k y — 0. Top panel 
shows the band structure for the weak topological insulator with the 
two Dirac points at the Brillouin zone boundaries. Band structure for 
the strong topological insulator with the Dirac point inside the band 
gap (bottom panel). 



We see that for the case of strong topological insulator there 
appears a Dirac point in the gap in the middle of the Brillouin 
zone (BZ). For the set of parameters giving a weak topological 
insulators, there are two Dirac points located at the edges of 
the BZ. We note that the Dirac node in the spectrum exists 
not only for the simple cubic unit cell, but also for the more 
complicated fee- and bec-unit cells. 



D. entanglement entropy and spectrum 



also three. A conventional band insulator has all of the four 
indices Iqti = I^ti = ^wti = ^wti = +1 or equivalently 
(0;0,0,0). An index I = (— 1) {v a = 1) indicates a Z 2 topo- 
logical state with the odd number of surface Dirac modes. In 
a KI the symmetry index 8 m of a particular high symmetry 
point m is negative provided £k« < £/ is lower the f-energy 
ef. Thus if £k^=o < £/ at the T point, while fk^o > £/ 
for all other symmetry points, then /gxi = — 1, and hence the 
Kondo insulating state is a strong-topological insulator, ro- 
bust against disorder Fig. [T] Weak-topological insulators and 
topologically trivial insulators can in principle be found for 
different band structures and different values of £/. A partic- 
ularly interesting possibility is to tune topological phase tran- 
sitions between different types of insulators (e.g., by applying 
a pressure). Although we have been specifically considering a 
tight-binding model with a primitive unit cell, all our conclu- 
sions apply directly to systems adiabatically connected to this 
model. 



In this Section we independently re-derive our results from 
the previous subsections by employing the concept of the en- 
tanglement entropy. In discussions on topological insulators 
without electron-electron interactions it is implicitly assumed 
that the presence of the gapless edge modes is a signature of 
the topologically non-trivial insulating state. In fact, this as- 
sumption is confirmed within our description of Kondo in- 
sulators. It is interesting, however, to check the topological 
properties of our model by discussing the properties of the 
eigenfunctions only. Such an approach has been pioneered by 
Freedman and collaborators 40 who showed that topologically 
nontrivial states of matter can exist without exhibiting the chi- 
ral edge modes. In this and the following Section we will dis- 
cuss in detail the topological properties of the eigenfunctions 



governed by our effective model Hamiltonian (21 



As it has been extensively discussed in the literature (for the 
more recent accounts see ET1 - I431 and references therein), en- 
tanglement entropy can be used to distinguish the topological 
phases from the non-topological ones. The following criterion 
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is used: the topologically nontrivial state should have non- 
zero entanglement entropy when the latter can not be tuned to 
zero by an adiabatic change of the parameters of the system 43 . 

As an aside, we note that our effective Hamiltonian ( pi) 
is a single-particle Hamiltonian and therefore, by calculating 
its entanglement spectrum we can also test the idea of adia- 
batic connectivity between our interaction-driven topological 
Kondo insulators and non-interacting topological insulators. 
The latter, however, cannot be adiabatically connected to triv- 
ial band insulators without making the system gapless. Note, 
for the trivial insulators we adopt the following definition 43 : 
upon adiabatic change in the hopping elements to zero, a triv- 
ial insulator goes into an atomic insulator without closing the 
energy gap along the adiabatic path. 

The entanglement entropy can be generally written as 

Sent = ~J2 l0 §^ + (! - &) - GO) ' < 24 > 

a 

where {£ a } are the single-particle entanglement eigenvalues, 
subscript a labels the eigenstates. When the entanglement 
eigenvalues are neither zero or one, the entanglement entropy 
is non-zero and therefore the system is topologically non- 
trivial. In particular, for translationally invariant topological 
insulator the spacial cut reveals the surface states and yields 
non zero entanglement entrop y 42 ! 43 1 i n t±iis case the entan- 
glement eigenvalues are also labeled by the conserved com- 
ponents of the momentum, say, £ a (k x ,k y ) for the cut in the 
xy-plane. Thus the problem of checking whether the insula- 
tor is topological or not reduces to the problem of determining 
the entanglement eigenvalues. For the case of the Kondo in- 
sulators, the computation of the entanglement spectrum may 
serve as an additional indicator of the nontrivial nature of 
their ground state especially for the case of complicated lat- 
tice structure when the simple approaches for the computation 
of the Z 2 indices do not apply. The procedure of how these 
eigenvalues are computed will be given below. 

In this Section we will evaluate the entanglement entropy 
for the mean field Hamiltonians pi) , using the Peschel's 
The entanglement spectrum is determined by 
correlation function 



G 



a/3 



(25) 



where ip ia creates an electron in state a = 1, .... 4 (conduc- 
tion or /- electron with spin up or down) on site i and the 
expectation value is evaluated in the ground state. Introduc- 
ing the normal operators 7„k, where n is the number of the 
eigenvalues: 



ipta = ^ elk riu na(k)7nk, 



(26) 



where Nf, = 2 is the number of the occupied bands and 
M na (k) are the eigenvectors. For the correlation function we 
find 

2 

G?f = £ e ik ^->) <a(kK/?(k) (27) 

k n=l 



summation goes over all components of the momentum k. 

Let us imagine now that our system is cut in two halves 
along a given spacial directions. To be specific, let us make 
the cut along the xy-plane, so that k x and k y are conserved. 



The entanglement spectrum £ a (kj_), kj^ = (k x , k y ) Eq. (24i 



will then be given by the eigenvalues of the following matrix: 



Gtf(k X ,ky)=j2 e ' 



k z -(z 



2 

2 ^<aW^W (28) 

n=l 



where i, j are confined to the right (or left) part of the sys- 
tem. Specifically, we need to solve the following eigenvalue 
problem: 



E 

3,P 



(kx) (29) 



We show the results of our computation for the model 
Hamiltonian pi) with Nb = 2 on Figs. [3] and [4] As we can 



see, depending on the position of renormalized /-level rela- 
tive to the bottom of the conduction band, we find either singe 
node or two nodes in the entanglement spectrum. The number 
of the nodes is equal to the number of nodes of the surface 
states inside the insulating gap, in complete agreement with 
our expectations. 

To summarize, our results from this Section confirm that 
Z 2 -odd topological Kondo insulators cannot be adiabatically 
connected to Z 2 -even insulators adiabatically without vanish- 
ing of the insulating gap along the adiabatic path. At the same 
time we see one-to-one correspondence between the Kondo 
insulators and non-interacting Z 2 topological insulators con- 
firming the idea of adiabatic connectivity between the two dis- 
cussed in the Introduction. 



IV. CONSTRUCTION OF THE WANNIER WAVE 
FUNCTIONS 

Within our model of Kondo insulators we can also ad- 
dress the problem of constructing maximally localized Wan- 
nier functions (WF). This question has several importa nt ap- 
plications in the general theory of topological insulator d 45 ! 46 ! 
such as calculation of the Z 2 indices as well as characteriza- 
tion of the topological structure using first principles calcula- 
tions. 

Generally, the construction of the WF proceeds in two 
stages. The first stage has to do with the initial choice of the 
basis set before specifying a particular choice of the gauge. 
This needs be done in order to make the WF nonsingular 
across the whole Brillouin zone. The gauge is then fixed by 
imposing the certain criterion. As an example, maximum lo- 
calization criterion is typically used 47 . Apart from the prob- 
lem related to the arbitrariness in the choice of Wannier func- 
tions, there exists a topological obstruction for constructing 
the Wannier functions for Chern insulators realized in systems 
with broken time-reversal symmetr y 48 ! 49 ! As it turns out, in 
the case of the Z 2 topological insulators there is also a topo- 
logical obstruction albeit a less severe one. As it was recently 
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FIG. 3: Entanglement spectrum for the model of Ce-based Kondo 
insulator: (a) weak topological insulator (0; 111), (b) strong topo- 
logical insulator (1; 000) and (c) trivial insulator (0; 000). For the 
presentation purposes we choose only two eigenvalues for each mo- 
mentum k i . 



discussed in Ref. [45 1 for the Kane-Mele model, it is impos- 
sible to construct the time-reversal invariant basis set of the 
Wannier functions, but one can construct the basis set consist- 
ing of the non-Kramers pairs. The above mentioned arbitrari- 
ness in the definition of the WF is then fixed by the criterion 
of maximum localization 4 -^! 

In this section we will specifically apply the prescription 
developed in Ref. P31 to construct the basis set which then 
can be used to initialize the procedure to compute the maxi- 
mum localized WF. On one hand, this should provide another 
example of the manifestation of above mentioned obstruction 
and the way it can be resolved. On the other hand, it gives an 
insight into the structure of the wave functions describing the 
quasiparticles in the occupied bands. 



A. preliminaries 

Below we will follow almost verbatim the discussion in 
Refs. B45I47I . For variety of applications (i.e. numerical 
calculations) it is required that the Bloch-like wave functions 
must remain smooth across the whole Brillouin zone (BZ). 
The problem is that the specific choice of the Bloch functions 
\ipnk) 1S not unique, since these wave functions have an ad- 
ditional gauge freedom originating from possibility of mixing 
with the wave functions describing the occupied bands: 



mk/ 



(30) 



(here the summation goes over the occupied bands). For all 
the practical purposes, however, the freedom of choosing the 
proper gauge transformation must be removed by applying 
some restrictions on choosing the specific gauge. The latter 
uses the criterion of maximum localization of the WF 4 -^. WFs 
are defined by 



W n (r - R) = 



n 



(2tt) 3 



ik R 



(31) 



BZ 



where Vt is a volume of the unit cell and V'nk(r) = (rlVVik) 
are the Bloch wave functions, n is a band index and R is a 
position of a lattice site. 



The unitary transformation ( 30 1 can be initialized using the 
following procedure. One first chooses the set of localized 
trial wave functions \t^) and then form a set of new basis 
functions 



\fjk) = iV'»k)(V ) nk|Tik)) 1=1, JV 



(32) 



where Af is the number of the occupied bands. Since this new 
basis set is not orthonormal, one can adopt a Lowdin proce- 
dure and form the overlap matrix 

S m „(k) = (f mk |f„ k ). (33) 

Now we can use Eqs. ( |32|33| l to form a set of Bloch-like states 



|V>nk>=£[^ 1/2 (k) 



(34) 



These states, albeit not eigenstates of the Hamiltonian, should 
be the smooth functions of the quasi-momentum k and are 
used to construct the localized set of the WFs: 



Wn(r-R) 



(2tt) 3 



-ik R 



VWO (35) 



BZ 



The above construction breaks of the determinant of the ma- 
trix £ mn (k) vanishes in some points of the BZ. Thus the 



problems consists in finding the proper set of trial states ( 32 1 
such that det[5(k)] ^ 0. Finally, we note that required de- 
gree of localization can be achieved by employing the iterative 
procedure^. 
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FIG. 4: Evolution of the entanglement spectrum with the change in the position of the renormalized /-level, s/. Top four panels show the 
change in the entanglement spectrum as the system goes from the trivial insulator to the strong topological insulator. Bottom four panels 
illustrate the changes in spectrum as system goes from strong to weak topological insulator. Note that when e / is exactly at the boundary 
separating different insulating phases so that the bulk insulating gap vanishes, the vertical lines in the entanglement entropy spectrum reflect 
the absence of the surface modes at the boundaries. 



We now construct the Wannier functions for our mean field 



model described by the Hamiltonian (21 1. For the Bloch wave 
functions we write 



IV'nk) = } C S nk\sk), 



s=l 



(36) 



where n = 1,2 labels the occupied bands, coefficients C snk 
are the components of the eigenvectors of the Hamiltonian 



(21 1 and the summation goes over the components of gen- 
eralized spinor which includes spinfull conduction (c) and f- 
elecron (f) states: 



| s = i,k) = 4 



3,k) 



fl t \o), 



= 2,k)=4,|0), 



4,k) 



kj.1 



(37) 



In Eq. (36 1 the basis functions |sk) are defined on the each 



site on the lattice R, i.e. 

|*k) = 



J2 e lk r S(r - R) 



(38) 



In what follows we adopt the method outlined above to con- 
struct the Wannier functions for our model Kondo insulators. 



B. choice of the basis 

Onset of the coherence in the Kondo lattice can be inter- 
preted as an emergence of new quasi-particles which are the 
linear superposition of the localized and conduction states. 
Since the newly formed quasiparticle band is narrow, the spec- 
tral weight is mostly governed by the /-states. Thus, to con- 
struct the Wannier functions we first consider the basis on /- 
states only: 



rik, 



= |3k), 



?2k 



|4k) 



(see Eqs. ( 37|38 1). For the new set of basis vectors (32i with 
the help of Eqs. ( 36|39 1 this implies 



Ink) = C'aiklV'ik) + C'lakl^k), 

|f 2k ) = C| lk |Vlk) + C 4 * 2k |^2k), 

For the determinant of the matrix 5(k) we find 



(40) 



det[S(k)] =(|C 31k | 2 + |C 32k | 2 )(|C 4 i k | 2 + |C 42k | 2 ) 

— |C 3 i k C| lk + C 32k C 42k | 2 



(41) 



We present the results on Fig. [5] We see that the determinant 
of the matrix ( 33 i is zero near the T-point which means that 
the choice ( 39 1 is not suitable for construction of non-singular 
Bloch functions and consequently Wannier functions. The 
same result holds for the trial basis built out of the conduc- 
tion states, | lk) and |2k) as well as their linear combinations. 

As we have mentioned above, formation of the coherence in 
the Kondo lattice can be seen as a formation of the new states 
(or quasiparticles) as a result of the hybridization between the 
conduction and / electrons. Motivated by this observation, let 
us try the following two trial basis wave-functions: 



|r lk ) = — (|lk) + |3k)), 
|r 2k ) = ~(|2k)-|4k)), 



(42) 



Note that the trial basis functions do not transform into each 
(39) other by time reversal operator, so they do not form a Kramers 
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FIG. 5: Plot of the dependence of det[S(k)] along the path in the 
BZ. The elements of the matrix S(k) has been obtained using the 
trial basis set, which consists of (a) non-Kramers pair of states each 
containing the superposition between the conduction and /-states; 
(b) Kramers pair of /-states and (c) Kramers pair states with linear 
superposition of conduction and /-electron wave functions. Deter- 
minant does not vanish anywhere in the BZ only for the basis (c). 



doublet. If follows 

|~ \ (C7lk + Qilk) 

Plk) = 7= 



lk 



|~ \ (Qilk ^41k) I j, \ 
= 71 



(^12k + ^32k) 
(^22k ~ ^42k) 

V2 



l^2k), 



l^2k), 



(43) 



The determinant of the matrix S(k) up to the numerical pre- 
factor is 



C 32k \ 2 )x 

2\ 



det[5(k)] =(|C uk + C 31k | 2 + |C 12k 

(|C21k — C4ik| 2 + |C22k — C*42k 
~ l(C*lk + CjllkX^lk - C 4 i k ) + 
+ (C*2k + C32k)(C22k - C 42 k)| 2 



(44) 



V. CONCLUSIONS 

In this paper we have discussed the conditions for the emer- 
gence of chiral surface states in semiconducting /-electron 
systems. We considered an insulating state in heavy fermion 
systems which appears at finite temperatures as a result of 
strong interaction between the conduction and the predomi- 
nantly localized /-electrons. Having started with the periodic 
Anderson lattice model, we considered the low-energy version 
of that model, which takes into account the effect of Hubbard 
repulsion between the /-electrons on the level of renormaliza- 
tions to the /-electron energy and hybridization amplitudes. 
The key ingredient of our model is momentum dependent hy- 
bridization amplitudes. The momentum dependence of the 
amplitudes originates from the strong spin-orbit coupling in- 
teraction on /-sites. The analysis of the topological structure 
of the newly formed insulating state is greatly simplified for 
the systems with simple cubic unit cell. In that case, the form 
factors vanish at high symmetry points of the BZ. This em- 
beds the topological singularities into the valence band, so that 
when the form-factors have p- or /-wave symmetry it imme- 
diately leads to the topological insulator. 

To describe the physics of Ce-based Kondo insulators, we 
considered the simplest model containing single conduction 
band hybridized with the Kramers doublet of /-states. We 
find that there will always be chiral surface states, when hy- 
bridization gap does not have nodes. The robustness of these 
states with respect to disorder is determined by the position 
of the renormalized /-level relative to the bottom of the con- 
duction band. We then verify our results for both models by 
calculating the entanglement entropy spectrum. Finally, we 
also discuss how to choose the basis for constructing Wan- 
nier wave functions, which are well defined everywhere in the 
Brillouin zone. It is interesting to note that the required ba- 
sis relies on superposition between the conduction and the lo- 
calized /-states. More importantly, this agrees with common 
view that a heavy quasiparticle is the quantum many-body su- 
perposition of conduction and /-states. 



We present the resulting dependence det[S'(k)] on momentum 
on Fig. [5] As we have expected, the determinant does not van- 
ish anywhere within the BZ which means we have succeeded 
in constructing the wave functions \ij) n k) (34i. In fact, we 



find det[5(k)] = 1 for the non-Kramers basis set (43 I. Fi- 
nally, we note that in agreement to the results of Ref. 14511 ob- 
tained for the Z2-odd phase in the Kane-Mele model, here our 
non-singular basis set also consists of the non-Kramers pair of 
states. To summarize, we have demonstrated that the basis for 
the Bloch wave functions can be chosen in such a way that no 
singularities are generated across the Brillouin zone. 
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